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Abstract. Recent analyses of helioseismological observa- 
tions seem to suggest the presence of two new phenomena 
connected with the dynamics of the solar convective zone. 
Firstly, there are present torsional oscillations with peri- 
ods of about 11 years, which penetrate significantly into 
the solar convection zone and secondly, oscillatory regimes 
exist near the base of the convection which are markedly 
different from those observed near the top, having either 
significantly reduced periods or being non-periodic. 

Recently spatiotemporal fragmentation / bifurcation 
has been proposed as a possible dynamical mechanism to 
account for such observed multi-mode behaviours in dif- 
ferent parts of the solar convection zone. Evidence for this 
scenario was produced in the context of an axisymmetric 
mean field dynamo model operating in a spherical shell, 
with a semi-open outer boundary condition and a zero or- 
der angular velocity obtained by the inversion of the MDI 
data, in which the only nonlinearity was the action of the 
Lorentz force of the dynamo generated magnetic field on 
the solar angular velocity. 

Here we make a detailed study of the robustness of this 
model with respect to plausible changes to its main ingre- 
dients, including changes to the a and 77 profiles as well 
as the inclusion of a nonlinear a quenching. We find that 
spatiotemporal fragmentation is present in this model for 
different choices of the rotation data and as the details of 
the model are varied. Taken together, these results give 
strong support to the idea that spatiotemporal fragmen- 
tation is likely to occur in general dynamo settings. 
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1. Introduction 

Recent analyses of the helioseismological data, from both 
the Michclson Doppler Imager (MDI) instrument on board 
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the SOHO spacecraft (Howe et al. 2000a) and the Global 
Oscillation Network Group (GONG) project (Antia & 
Basu 2000) have provided strong evidence to indicate that 
the previously observed time variation of the differential 
rotation on the solar surface - the so called 'torsional os- 
cillations' with periods of about 11 years (e.g. Howard 
& LaBonte 1980; Snodgrass, Howard & Webster 1985; 
Kosovichev & Schou 1997; Schou et al. 1998) - penetrates 
into the convection zone, to a depth of at least 9 percent 
of the solar radius. Torsional oscillations are thought to 
be a consequence of the nonlinear interactions between 
the magnetic fields and the solar differential rotation. A 
number of attempts have been made to model these os- 
cillations. These include modelling the dynamical feed- 
back of the large scale magnetic field on the turbulent 
Reynolds stresses that generate the differential rotation 
in the convection zone (the 'nonlinear A-effect': Kitchati- 
nov 1988; Riidiger & Kitchatinov 1990; Kitchatinov et 
al. 1994; Kiiker et al. 1996; Kitchatinov & Pipin 1998; 
Kitchatinov et al. 1999), as well as models in which the 
nonlinearity is through the direct action of the azimuthal 
component of the Lorentz force of the dynamo generated 
magnetic field on the solar angular velocity (e.g. Schuessler 
1981; Yoshimura 1981; Brandenburg & Tuominen 1988; 
Jennings 1993; Covas et al. 2000a,b; Durney 2000). 

Further studies of these data have produced intrigu- 
ing, but rather contradictory results. Howe et al. (2000b) 
find evidence for the presence of such oscillations around 
the tachocline near the bottom of the convection zone, but 
with markedly shorter periods of about 1.3 years. On the 
other hand, Antia & Basu (2000) do not find such oscilla- 
tions at the bottom of the convective zone. Whatever the 
true dynamical behaviour at these lower levels may be, the 
crucial point about these recent results is that they both 
seem to indicate the possibility that the variations in the 
differential rotation can have different dynamical modes of 
behaviour at different depths in the solar convection zone: 
oscillations with a very different period in the former case 
and non-periodic behaviour in the latter. 

Clearly, further observations are required to clarify this 
situation. However, whatever the outcome of such obser- 
vations, it is of interest to ask whether such different varia- 
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tions can in principle occur in different spatial locations in 
the convection zone and, if so, what could be the possible 
mechanism(s) for their production. This is of particular 
interest given the inevitable errors in hclioscismological 
inversions, especially as depth increases. 

Recently, spatiotemporal fragmentation/bifurcation 
has been proposed as a possible dynamical mechanism 
to account for the observed multi-mode behaviour at 
different parts of the solar convection zone (Covas et al. 
2000b) (hereafter CTM). This occurs when dynamical 
regimes which possess different temporal behaviours 
coexist at different spatial locations, at given values of 
the control parameters of the system. The crucial point 
is that these different dynamical modes of behaviour can 
occur without requiring changes in the parameters of the 
model, in contrast to the usual temporal bifurcations 
which result in identical temporal behaviour at each 
spatial point, and which occur subsequent to changes 
in the model parameters. Also, as we shall see below, 
spatiotemporal fragmentation/bifurcation is a dynamical 
mechanism, the occurrence of which does not depend 
upon the detailed physics at different spatial locations. 

Evidence for the occurrence of this mechanism was 
produced in CTM in the context of a two dimensional ax- 
isymmetric mean field dynamo model operating in a spher- 
ical shell, with a semi-open outer boundary condition, in 
which the only nonlinearity is the action of the azimuthal 
component of the Lorentz force of the dynamo generated 
magnetic field on the solar angular velocity. The zero order 
angular velocity was chosen to be consistent with the most 
recent helioseismological (MDI) data. Despite the success 
of this model in producing spatiotemporal fragmentation, 
a number of important questions remain. Firstly, there are 
error bars due to the nature of the observational data as 
well as inversion schemes used. Secondly, the model used 
by CTM is approximate and includes many simplifying 
assumptions. As a result, to be certain that spatiotempo- 
ral fragmentation can in fact be produced as a result of 
such nonlinear interactions, independently of the details 
of the model employed, it is necessary that it is robust, 
i.e. that it can produce such behaviour independently of 
these details. 

The aims of this paper are twofold. Firstly, we make a 
systematic study of of this mechanism by making a com- 
parison of the cases where the zero order rotational ve- 
locity are given by the inversions of the MDI and GONG 
data respectively. Given the detailed differences between 
these rotation profiles, this amounts to studying the ro- 
bustness of the mechanism with respect to small changes 
in the zero order rotation profile. 

Secondly, we study the robustness of this model (in 
producing spatiotemporal fragmentation) with respect to 
a number of plausible changes to its main ingredients, and 
demonstrate that the occurrence of this mechanism is not 
dependent upon the details of our model. 



We show that in addition to producing butterfly di- 
agrams which are in qualitative agreement with the ob- 
servations as well as displaying torsional oscillations that 
penetrate into the convection zone, as recently observed by 
Howe et al. (2000a) and Antia & Basu (2000), and studied 
by Covas et al. (2000a), the model can produce qualita- 
tively different forms of spatiotemporal fragmentary be- 
haviours, which could in principle account for either of 
the contradictory types of dynamical behaviour observed 
by Howe et al. (2000b) and Antia & Basu (2000) at the 
bottom of the convection zone. 

The structure of the paper is as follows. In the next 
section we outline our model. Section 3 contains our de- 
tailed results for both MDI and GONG data. In section 4 
we study the robustness of spatiotemporal fragmentation 
with respect to various changes in the details of our model 
and finally section 5 contains our conclusions. 

2. The model 

We shall assume that the gross features of the large scale 
solar magnetic field can be described by a mean field dy- 
namo model, with the standard equation 



dB 

~dt 



Vx(uxB + aB-t)VxB), 



(1) 



where B and u are the mean magnetic field and the mean 
velocity respectively. The quantities a (the a effect) and 
the turbulent magnetic diffusivity rj t , appear in the process 
of the parameterisation of the second order correlations 
between the velocity and magnetic field fluctuations (u' 
and B'). Here u = v<p— the term proportional to V77 
represents the effects of turbulent diamagnetism, and the 
velocity field is taken to be of the form v = vq + v' , where 
vo — f2o?~sin#, Qq is a prescribed underlying rotation law 
and the component v' satisfies 



cV (VxB)xB- _ 2 , 
at fioprsmt) 



(2) 



where D 2 is the operator |^ + f ^ + ^r^g ( W ( sin 9 W ) ~ 
-^s) and fi is the induction constant. The assumption 
of axisymmctry allows the field B to be split simply into 
toroidal and poloidal parts, B = By+Bp = V x A<p, 
and Eq. (||) then yields two scalar equations for A and B. 
Nondimensionalizing in terms of the solar radius R and 
time R 2 /r]o, where 770 is the maximum value of ij, and 
putting £1 = a = ao<5, i] = r/r>rj, B — £?oB and 

v' = Q*Rv' , results in a system of equations for A, B 

dA 

— = W r Bg 

or 
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where w = — The dynamo parameters are R a = 

a R/r)0i = ^*R 2 /va, p r = vo/m, and V = v/vo, where 
J7* is the solar surface equatorial angular velocity. Here vq 
and r)o are the values of the turbulent magnetic diffusivity 
and viscosity respectively and P r is the turbulent Prandtl 
number. The density p is assumed to be uniform. 

For our inner boundary conditions we chose B = 
on r — ro (which enforces angular momentum conserva- 
tion in the dynamo region, as discussed in Moss & Brooke 
2000). At the outer boundary, we used an open bound- 
ary condition dB/dr = on B, and vacuum boundary 
conditions for Bp (as in CTM). The motivation for the 
former condition is that the surface boundary condition 
is ill-defined, and there is some evidence that the more 
usual B = condition may be inadequate. This issue has 
recently been discussed at length by Kitchatinov, Mazur 
& Jardine (2000), who derive 'non- vacuum' boundary con- 
ditions on both B and Bp. 

Equations @, @ and (||) were solved using the code 
described in Moss & Brooke (2000) (see also Covas et al. 
2000a, b) together with the above boundary conditions, 
over the range J"o<r<l,O<0<7r. We set ro = 0.64; 
with the solar convection zone proper being thought to 
occupy the region r > 0.7, the region ro < r < 0.7 can 
be thought of as an overshoot region/tachocline. In the 
following simulations we used a mesh resolution of 61 x 
101 points, uniformly distributed in radius and latitude 
respectively. 

The model employed in CTM had the following ingre- 
dients. In the interval 0.64 < r < 1, fio was given by the 
inversion of the MDI data obtained from 1996 to 1999 by 
Howe et al. (2000a). Here, in addition to this, we shall for 
comparison also use the S7 given by the inversion of the 
GONG data obtained from 1995 to 1999 by Howe et al. 
(2000b). The form of a was taken as 



a = a r {r)f{6), 



(5) 



where f(8) = sin 4 9 cos 9 (cf. Riidiger & Brandenburg 
1995). The model possessed radial dependence in both a 
and the turbulent diffusion coefficient 77: the a profile was 
chosen by setting a r = 1 for 0.7 < r < 0.8 with cubic 
interpolation to zero at r = ro and r = 1, with the con- 
vention that a r > and R a < 0. The rj profile was chosen 
in order to take into account the likely decrease of rj in 
the overshoot region, by allowing a simple linear decrease 
in 77 from 77 = 1 at r = 0.8 to fj — 0.5 in r < 0.7. 

In the following sections we shall, in addition to these 
particular profiles, also consider variations to them in or- 
der to study the robustness of the model (in produc- 
ing spatiotemporal fragmentation) with respect to such 
changes. Also for the sake of comparison, we shall use the 
interpolation on the GONG data for the rotation law as 
well as allowing other changes. 



3. Torsional oscillations and fragmentation using 
the MDI and GONG data sets 

In this section, we study the presence of spatiotemporal 
fragmentation and torsional oscillations, using the zero 
order angular velocity obtained by the inversion of the 
GONG data (see Fig. |). This allows a detailed compari- 
son to be made with the results of CTM which employed 
the corresponding MDI data (see Fig. ^), including the 
magnetic field structure and strengths as well as the na- 
ture of torsional oscillations, as a function of depth and 
latitude. 




Fig. 1. Isolines of the time average of the angular veloc- 
ity of the solar rotation, obtained by inversion techniques 
using the GONG data (Howe et al. 2000b). Contours are 
labelled in units of nHz. 




Fig. 2. Isolines of the time average of the angular veloc- 
ity of the solar rotation, obtained by inversion techniques 
from the MDI data (Howe et al. 2000a). Contours are la- 
belled in units of 11Hz. 
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To begin with, we calibrated our model in each case so 
that near marginal excitation the cycle period was about 
22 years. This determined i? w to be 44000, corresponding 
to ?7o ~ 3.4 x 10 11 cm 2 sec -1 , given the known values of 
VL* and R. 

In both cases, the first solutions to be excited in the 
linear theory are limit cycles with odd (dipolar) parity 
with respect to the equator, with marginal dynamo num- 
bers R a « —3.05 and —3.11 respectively for the GONG 
and MDI data. The even parity (quadrupolar) solutions 
are also excited at similar marginal dynamo numbers of 
R a —3.20 and —3.22, respectively. We note that these 
marginal dynamo numbers are slightly different from those 
reported in CTM as we have used here f(9) = sin 4 9 cos 9 
for the ^-dependence of the a effect, whilst CTM used the 
prescription f(9) = sin 2 9 cos 9. 

Also, to be consistent with CTM, we chose the value 
of the Prandtl number in this section to be P r = 1.0. For 
the parameter range that we investigated, the even parity 
solutions can be nonlinearly stable. Given that the Sun 
is observed to be close to an odd (dipolar) parity state, 
and that previous experience shows that small changes in 
the physical model can cause a change between odd and 
even parities in the stable nonlinear solution, we chose to 
impose dipolar parity on our solutions, effectively solving 
the equations in one quadrant, < 9 < 7r/2, and imposing 
appropriate boundary conditions at the equator 9 = n/2. 

With these parameter values, we found that this 
model, with the underlying zero order angular velocity 
chosen to be consistent with either the MDI or the GONG 
data, is capable of producing butterfly diagrams which are 
in qualitative agreement with the observations, as can be 
seen in Figs. || and [|. We note that the butterfly diagrams 
in the GONG case tend to be smoother and to have weaker 
polar branches, presumably because the MDI inversion is 
respectively less smooth and has a distinctive polar jet 
close to r = 0.95i? Q . 

The model can also produce torsional oscillations in 
both cases (see Figs. || and ||), that penetrate into the 
convection zone, in a manner similar to those deduced 
from recent helioseismological data (Howe et al. 2000a; 
Antia & Basu 2000) and studied in Covas et al. (2000a, b). 

We also found the model to be capable of producing 
spatiotemporal fragmentation near the base of the convec- 
tion zone, i.e. we found there different dynamical modes 
of behaviour in the differential rotation, including oscilla- 
tions with reduced periods, as well as non-periodic vari- 
ations as in CTM. These coexist with the near-surface 
torsional oscillations of the form shown in Figs. || and ^). 
To demonstrate this, we have plotted in Figs. and || the 
radial contours of the angular velocity residuals 5fl as a 
function of time for a cut at latitude 30° . To demonstrate 
the period halving produced by the fragmentation more 
clearly, we have also plotted in Fig. || the residuals 6fl at 
different values of R a . 




Time (years) 



Fig. 3. Butterfly diagram of the toroidal component of 
the magnetic field B at fractional radius r = 0.95 for the 
MDI data. Dark and light shades correspond to positive 
and negative values of B<p respectively. Parameter values 
are R a = -11.0, P r = 1.0 and R u = 44000. 
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Fig. 4. Butterfly diagram of the toroidal component of 
the magnetic field B at fractional radius r = 0.95 for the 
GONG data. Dark and light shades correspond to positive 
and negative values of B^ respectively. Parameter values 
are R a = -9.0, P r = 1.0 and R u = 44000. 

We also find that, in all cases, torsional oscillations 
with the same period persist in all the regions above the 
fragmentation level as can, for example, be seen from Figs. 

and ||. The butterfly diagrams for the toroidal magnetic 
field, on the other hand, keep the same period and quali- 
tative form throughout the convection zone, including the 
spatial regions that experience fragmentation, as shown 
in Figs, [h] and [ll]. Thus the fragmentation in the angular 
velocity residuals do not seem to be present in the mag- 
netic field. Fig. [ll] shows clearly that the fragmentation 
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Fig. 5. Variation of the perturbation to the zero order 
rotation rate in latitude and time, revealing the migrating 
banded zonal flows, taken at fractional radius r = 0.95, 
with the MDI data. Parameter values are as in Fig. || 




5 10 

Time (years) 



20 



Fig. 6. Variation of the perturbation to the zero order 
rotation rate in latitude and time, revealing the migrating 
banded zonal flows, taken at fractional radius r = 0.95 
with the GONG data. Parameter values are as in Fig. |[ 



only occurs in 5fl and that the magnetic field retains its 
typical 22 year cycle throughout the dynamo region. 

We also made a detailed study of the magnetic field 
evolution and the magnitude of torsional oscillations in 
each case. Fig. |l2] summarises the results of the average 
magnetic energy as a function of R ai for both the MDI 
and the GONG data, and demonstrates that both data 
sets produce almost the same average magnetic energy 
for a given R a . Similarly, we have plotted in Fig. 13, the 



maxima (of the absolute value) of the residuals of the dif- 
ferential rotation. Sfl, as a function of R a , for both MDI 



and the GONG data, which shows that for large values of 
|i? a |, the GONG data produce larger residuals. 

We end this section by summarising the qualitative 
modes of spatiotemporal behaviour our detailed numeri- 
cal results have produced, for both MDI and GONG data. 
Our results show that our model is capable of producing 
three qualitatively different spatiotemporal modes of be- 
haviour: (i) regimes in which there is no deformation of the 
torsional oscillation bands in the (r, t) plane, and hence no 
changes in phase or period through the convection zone; 
(ii) regimes in which there is deformation in the oscillatory 
bands in the (r, t) plane, resulting in changes in the phase 
of the oscillations, but no changes in their period and (hi) 
regimes with spatiotemporal fragmentation, resulting in 
changes both in phase as well as in period/behaviour of 
oscillations. An example of the difference between such 
regions is given in Fig. 0, which shows a spatiotempo- 
ral fragmentation resulting in period halving, as well as 
a phase change between the oscillations near the surface 
and those deeper down in the convection zone. We should 
emphasise that by period we mean here the time between 
the maxima (minima) , rather than time between repeated 
sequences. We note that given the presence of noise, the 
error bars induced by the inversion and the shortness of 
the observational interval, the time between the maxima 
(minima) is the most relevant quantity from the point of 
view of comparison with observations. 

Our results also show that using the zero order rotation 
produced by both the MDI and GONG data gives rise to 
qualitatively similar results regarding both the spatiotem- 
poral fragmentation and the torsional oscillations. There 
are, however, detailed quantitative differences. An exam- 
ple of this relates to the case (ii) above, where the location 
of where the fragmentation occurs is higher for MDI than 
for GONG at higher values of \R a \ (see Fig. |l4|). Further- 
more, the phase shift is less pronounced for the MDI case 
than for the GONG case, specially at higher |i? a | values 
(see Fig. |l|). 

4. Robustness of the model 

As noted above, the model of the solar dynamo used by 
CTM to demonstrate the presence of torsional oscillations 
and spatiotemporal fragmentation is inevitably approxi- 
mate in nature and contains major simplifications and pa- 
rameterizations. It is therefore important to ask to what 
extent the spatiotemporal fragmentation, as well as the 
qualitative features of torsional oscillations found in CTM, 
depend on the details of the model employed. 

In order to go some way towards answering these ques- 
tions, we shall in this section study the robustness of this 
model with respect to plausible changes to its main ingre- 
dients, considering these changes one by one. 

A feature of the model considered by CTM is the pres- 
ence of radial variations in both a and 77, and in particular 
the fact that a was taken to be zero in r < 0.70. An obvi- 
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Fig. 7. Radial contours of the angular velocity residuals 
Sft as a function of time for a cut at latitude 30° , with the 
MDI data. Parameter values are as in Fig. ||. 
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Fig. 8. Radial contours of the angular velocity residuals 
5ft as a function of time for a cut at latitude 30°, for the 
GONG data. Parameter values are as in Fig. 0. 



ous question is whether this variation is the source of the 
fragmentation described above, for example does it cause 
the dynamo region essentially to consist of two disparate 
parts, each with its own set of properties? This is per- 
haps unlikely, as there is no independent dynamo action 
in r < 0.7, as a = there. Nevertheless, to resolve this 
issue, we study the effects of changes in the radial varia- 
tions of both a and rj and subsequently examine a model 
in which neither vary radially. 

Given the qualitative similarity of the results produced 
using the MDI and the GONG data, and in order to keep 
the extent of the computations within tractable bounds, 
we consider only the MDI data in the following sections. 



1.0 : 




10 20 30 
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Fig. 9. 'Period halving' at r = 0.68 and latitude 30°. The 
panels correspond, from top to bottom, to R a = —6.0 and 
— 11.0 respectively, and display increasing relative ampli- 
tudes of the secondary oscillations. The remaining param- 
eter values are as in Fig. [|. 
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Fig. 10. Butterfly diagram of the toroidal component of 
the magnetic field B at fractional radius r = 0.68 for the 
MDI data, showing no signs of fragmentation. Parameter 
values are as in Fig. 

Robustness with respect to changes in the a profile 

We studied this by setting a r = 1 throughout the compu- 
tational region, which changes the a profile substantially. 
In spite of this, we found qualitatively similar forms of 
butterfly diagrams, torsional oscillations and spatiotem- 
poral fragmentation in this case as in the case of CTM. 
An example of such fragmentation in the radial contours 
of SQ is given in Fig. fl6[ 

We note that the change in the r-dependence of a does 
not change the calibration of R u required to produce a cy- 
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Fig. 11. Radial contours of the toroidal component of the 
magnetic field B as a function of time for a cut at latitude 
30°, for the MDI data, showing no signs of fragmentation. 
This is the counterpart of the Fig. || and clearly shows that 
the fragmentation occurs only for 6Q(r, 9, t) and that the 
magnetic field retains its typical 22 year cycle. Parameter 
values are as in Fig. |[ 



u 
W 

a 

5o 



T3 

CD 

u 
> 

< 



16 



14 



12 



10 




-3 -5 



-7 



MDI 
GONG 

-11 -13 -15 



R 



Fig. 12. The average magnetic energy 2?m as a function 
of i? Q , for both MDI and the GONG data. The parameter 
values are P r — 1.0 and iL, = 44000. 
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Fig. 13. Comparison of the maximum of 5fl(r = 
0.95i? Q ,6>) for both the MDI and the GONG data sets. 
The parameter values are P r = 1.0 and R u = 44000. 




Fig. 14. Variation of fragmentation level (the maximum 
radius for which more than one fundamental period exists) 
as a function of the control parameter R a for both the 
MDI and the GONG data sets. The parameter values are 
P r = 1.0 and = 44000. 



cle period of 22 years. It does, however, seem to change the 
details of the surface torsional oscillations, making them 
somewhat less realistic. It also, importantly from the point 
of view of observations, enhances the phase shift between 
the torsional oscillations at the top and below the defor- 
mation, to a phase change of <j> = —tt as can be seen from 
Fig. g7|. Note, however, that the phase shift <j> changes 
continuously with R a . 



4-2. Robustness with respect to r\ profile 

To obtain some idea of the sensitivity of the model to 
the chosen radial dependence of the turbulent diffusion 
coefficient rj, we then examined a model in which i){r) = 
1 throughout. In this way we removed from the model 
the inhomogeneity associated with rj at the base of the 
convection zone. 

We found that with this form of 77, R u needs to be 
recalibrated to R^ = 60000 in order to produce toroidal 
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Fig. 15. Variation of phase shift (ft as a function of 
the control parameter R a for both the MDI and the 
GONG data sets. The phase shift evaluated by compar- 
ing the residual time series Sfl(r = O.95i?0, 6 = 30°) and 
SQ(r — 0.75Rq,0 — 30°). A positive phase shift (ft in- 
dicates a forward shift of the torsional oscillation bands 
(top to bottom) and a negative phase shift an opposite 
one. The parameter values are P r = 1.0 and i? w = 44000. 
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Fig. 16. Radial contours of the angular velocity residuals 
Sfl as a function of time for a cut at latitude 30°. The 
parameter values are R^ = 44000, R a = —11.0 and Pr = 
1.0 with a r (r) = 1. 



field butterfly diagrams at the surface with a period of 22 
years. 

With this value of R^, we again found qualitatively 
similar forms of spatiotemporal fragmentation in this case, 
to those obtained by CTM. An example is given in Fig. 




We note that, compared with the behaviour of the 
original model depicted in Figs. ||, this model produces 
more realistic torsional oscillations at the surface than the 
model with the radial dependence of a removed. However, 
the butterfly diagrams of the toroidal magnetic field near 
the surface look slightly less realistic than those for the 
original model shown in Fig. |^. But, it must be borne in 
mind that in the absence of a definitive sunspot model, it 
is not obvious at which depth the toroidal field patterns 
correspond to the observed butterfly diagram. Further- 
more, as a result of the removal of the radial variation of 
r\ in this model, the effective magnetic turbulent diffusiv- 
ity changes and consequently the range of R a values for 
spatiotemporal fragmentation are somewhat higher than 
for the model employed by CTM. We note also that the 
fragmentation is now present at higher radius and is in 
this case detached from the lower boundary. 




Time (years) 



Fig. 17. Radial contours of the angular velocity residuals 
Sil as a function of time for a cut at latitude 30°. The 
parameter values are R^ = 60000, R a — —16.0 and Pr = 
1.0 with fj(r) = 1. 



4-3. Robustness to simultaneous removal of radial 
dependence of a and r\ 

In the last two sections we studied the robustness of the 
spatiotemporal fragmentation and the torsional oscilla- 
tions with respect to the removal of the radial dependence 
of a and n in turn. Interestingly, in both cases we found 
the spatiotemporal fragmentation to persist. Here we in- 
vestigate a model in which both a and n are independent 
of radius, so that (unrealistically) the convective region 
and overshoot layer are only distinguished by the nature 
of the given zero order rotation law. Thus we considered 
a modification of the model considered by CTM, given by 
taking a r = fj(r) = 1. Perhaps surprisingly, we still found 
spatiotemporal fragmentation, as can be seen in Fig. [l8|. 
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We note, however, that in this case both the butterfly 
diagrams of the toroidal field and the torsional oscillations 
near the surface look less realistic than the original model 
of Sect. |^. Our point is however that spatiotemporal frag- 
mentation appears to be a general property of the type 
of dynamo models studied, and is not crucially dependent 
on imposed spatial structures in the coefficients a and rj, 
that might be thought to distinguish the lower part of the 
dynamo region from the subsurface layers. We also note 
that again the fragmentation region is detached from the 
lower boundary (suggesting that this is a feature related 
to the rj profile) and extends to an even greater height, 
r = O.71i?0. It is interesting that, within our limited ex- 
perimentation, the details of the fragmentation, but not 
its existence, seem to be quite sensitive to the rj profile. 

We also studied the average magnetic energy and the 
residuals, SQ, as a function of R a with different a and rj 
profiles and these are depicted in Figs, [lj] and As can 
be seen, the values of <50 are larger for the case with radial 
variations of both a and rj present, whereas the average 
magnetic energy is largest for the case with uniform rj. 




Time (years) 



Fig. 18. Radial contours of the angular velocity residuals 
5Vt as a function of time for a cut at latitude 30°. Param- 
eter values are i? w = 60000, R a = —15.0 and Pr = 1.0 
with r)(r, 9) = 1 and a r (r) = 1. 
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Fig. 19. Comparison of the average magnetic energy Em 
as a function of R a for models with different a profiles. 
The parameter values are P r = 1.0 and = 44000. The 
continuous curve represents the results of the model with 
a r = f(r), r)(r, 9) — r](r), the long-dashed curve the model 
with a r = 1 = n(r,9), the short-dashed curve the model 
with rj(r, 9) = 1 and the dotted-dashed curve the model 
with a r = 1. 




Fig. 20. Comparison of the maximum of SQ(r = 
0.95Rq,9) for models with different a profiles. The pa- 
rameter values and the characterisation of the curves is as 
in Fig. |TJ. 



4-4- Robustness with respect to a quenching 

Even though the model used by CTM was nonlinear (via 
the Lorentz force and the subsequent v' term in the in- 
duction equation), the magnitude of a was fixed ab initio. 
We note that the form of a r in the original model was 
initially chosen to represent implicit strong ev-quenching 
occurring in the overshoot layer, in that a = in r < 0.7. 



In this section we study the effects of having an additional 
nonlincarity in the form of an a quenching given by 



a r (r) 



sm 



' COSf 



1 + .9IB 2 



(6) 



where g is a quenching factor, in addition to the nonlin- 
carity given by Eq. <M). 
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In this connection it may be noted that there is an 
ongoing controversy regarding the nature and strength of 
a-quenching. This is not the place to rehearse the argu- 
ments for and against 'strong' alpha-quenching; the issue 
is unresolved, and we use (|^) as a commonly adopted non- 
linearity. 

We found that for values of g < 0.1 (depending some- 
what on R a and P r ), the model continues to produce tor- 
sional oscillations and spatiotemporal fragmentation. 

4-5. Robustness with respect to the Prandtl number P r 

In the results obtained by CTM, the Prandtl number was 
taken as P r = 1. In this section we study the effects of 
changing P r on the nature of spatiotemporal fragmenta- 
tion and torsional oscillations. 

Our studies show that, for our model, both spatiotem- 
poral fragmentation and torsional oscillations persist for 
values of Prandtl number given by P r > 0.4 (depending 
on R a ). Around P r ~ 0.4, however, a sudden transition 
seems to occur, such that below this value spatiotemporal 
fragmentation as well as coherent surface torsional oscil- 
lations seem to be absent. 

We also made a detailed study of the magnetic field 
strength as well as of the magnitude of the torsional oscil- 
lations, for different values of Prandtl number P r , Figs, ^l] 
and |22| show the behaviours of the average magnetic en- 
ergy and the residuals of the differential rotation, Sfl, for 
fixed for different values of Prandtl number P r . This 
is in agreement with the results of Kuker et al. 1996, who 
studied solar torsional oscillations in a model with mag- 
netic quenching of the Reynolds stresses (the 'A-effect'), 
with both studies showing a linear growth of the ampli- 
tude of the torsional oscillations as a function of (small) 
Prandtl number. In their case the amplitudes of the tor- 
sional oscillations for P r < 1 are weaker than observed, 
whereas in our model they are within physically more rea- 
sonable ranges, Sil « 0.5-1 nHz. However, the restriction 
to uniform density makes it uncertain how important these 
amplitude differences really are. 

Finally we also studied the average magnetic energy 
and the maximum value of the residuals, <50, as a func- 
tion of R a for a number of cases, including a nonlinearly 
quenched a, as shown in Figs. ^3] and ^4|. The decrease in 
P r results in the reduction of the average magnetic energy 
and the residuals, as can be seen from these figures. 

4-6. Variations in fragmentation level and phase 

Finally, in this section we briefly summarise our detailed 
results concerning the variations in the spatiotemporal 
fragmentation level (i.e. the maximum radius at which 
more than one period exists) as well as changes in the 
phase <j), as a function of the control parameters R a and 
P r for models with different a and rj profiles as well as 
those with nonlinear a-quenching. 




Prandlt number P r 

Fig. 21. The average magnetic energy Em as a function 
of P r . The parameter values are i? w = 44000, R a — —15.0 
and Pr = 0.5. 
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Fig. 22. The maximum of 8Q(r = 0.95i?o, 9) as a func- 
tion of P r when R u = 44000, R a = -15.0. 

Fig. |25| shows the radius of the upper boundary of the 
region of spatiotemporal fragmentation as a function of R a 
for different a and n profiles at P r = 1. As can be seen, 
for intermediate values of \R a \, the models with uniform 
n have a higher fragmentation level. Also they also have 
fragmentation at values of \R a \ significantly larger than 
for the model with a radial dependence of r\. The higher 
upper boundary of the fragmentation region for such mod- 
els is associated with the position of the fragmented cells, 
which are now detached from the bottom boundary, un- 
like the basic models described earlier (cf. Figs. 0, [l?] and 
|l8| ). At small values of \R a \ the removal of the radial de- 
pendence on the a profile does not change the maximum 
height at which fragmentation occurs, and even at higher 
\R a \ it does so only slightly. 

We have also studied the position of the boundary of 
the fragmentation region in models with cv-quenching as 
well as lower values of P r and the results are shown in Fig. 
. It can clearly be seen that the quenching of a changes 
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Fig. 23. Comparison of the average magnetic energy Em 
for models with a-quenching and smaller P r . In this fig- 
ure, R u = 44000, the continuous curve corresponds to 
the model with a r = f(r),r)(r,6) — r](r),P r = 1.0, the 
dot-dashed line is for the model with a-quenching and 
the dashed curve for the model with P r = 0.5 and no 
a-quenching. 
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Fig. 24. Comparison of the maximum of 50(r = 
O.95i?0,6*) for models with a-qucnching and smaller val- 
ues of P r . In this figure, R u = 44000, the continuous curve 
corresponds the model with a r = f(r),rf(r,6) = i](r), 
P r = 1.0, the dot-dashed line represents the model with 
a-quenching and the dashed curve the model with P r = 
0.5. 



significantly the position and amplitude of the fragmen- 
tation. Not only does it start at higher |i? a | but it also 
extends closer to the bottom boundary. The dependence 
on the Prandtl number is somehow less pronounced. Frag- 
mentation again sets in starts at higher values of \R a \ but 
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Fig. 25. Variation of spatiotemporal fragmentation level 
for different a and ij profiles. Here Cl(r, 0) is given by the 
MDI data, P r = 1.0, R u = 44000 and the continuous curve 
gives the results for the model with a r = f(r),r](r,6) = 
7](r), the long-dashed curve the model with a r = r)(r, 6) = 
1, the short-dashed curve the model with rj(r,0) = 1 and 
the dotted-dashed curve the model with a r = 1 (slightly 
displaced downwards for clarity). 
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Fig. 26. Variation of fragmentation level as a function of 
R ai in the presence of a-quenching as P r is reduced. Here 
Ra, = 44000 and the continuous curve gives the results 
for the model with a r = f(r),r](r,0) = r)(r),P r = 1.0 
(slightly displaced upwards for clarity), the dot-dashed 
line represents the model with a-quenching and the 
dashed curve the model with P r = 0.5. 
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the position of the fragmentation region is almost the same 
as that of the basic model with P r = 1. 

10 i 




Fig. 27. Variation of the phase shift as a function of R a 
for different a and r\ profiles. Here f2(r, 9) is given by the 
MDI data, P r — 1.0, R u — 44000 and the continuous curve 
gives the results for the model with a r = f(r),rj(r,9) = 
rj(r), the short-dashed curve the model with r](r,9) = 1 
and the dotted-dashed curve the model with a r — 1. We 
were unable to calculate uniquely the phase shift for the 
case without radial dependence in both the a and r\ pro- 
files and so we have omitted it. 

Further, we studied the shift in phase, 4>, between the 
oscillations near the top and bottom of the dynamo region. 
Fig. ^7] shows the phase shift as a function of R a , for 
different a and r\ profiles. As can be seen, the phase shift 
can be negative or positive. For the models with a r = 1, 
4> can be either positive or negative, while for the other 
models it is always positive. We note that the phases shifts 
for small enough \R a \ are correspondingly small, that is 
the torsional oscillations are all in phase throughout the 
convection zone. 

In Fig. |28| we have plotted the phase shifts for the a- 
quenched models. This shows that a-quenching does not 
dramatically change the behaviour (and the sign) of the 
phase shifts (apart from reducing it somewhat), for nearly 
all values of R a . 

5. Discussion 

We begin by acknowledging the considerable uncertain- 
ties and simplifications which are associated with all solar 
dynamo models, including our own. In particular, in the 
present context our assumption of uniform density may 
be important. 

We have made a detailed study of a recent proposal 
in which the recent results of helioseismic observations re- 
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Fig. 28. Variation of phase shift as a function of R a 
in the presence of a-quenching. Here R u — 44000, the 
solid line gives the results for the model with a quenching 
and the dashed line those for the model with a-quenching 
with g = 10~ 2 . 

garding the dynamical modes of behaviour in the solar 
convection zone are accounted for in terms of spatiotem- 
poral fragmentation. Originally, support for this scenario 
came from the study of a particular two dimensional ax- 
isymmetric mean field dynamo model operating in a spher- 
ical shell, with a semi-open outer boundary condition, 
which inevitably involved a number of simplifying and 
somewhat arbitrary assumptions. To demonstrate, to a 
limited extent at least, the independence of our proposed 
mechanism from the details of our model, we have shown 
that this scenario is robust with respect to a number of 
of plausible changes to the main ingredients of the model, 
including the a and rj profiles as well as the inclusion of 
nonlinear quenching and changes in the Prandtl number. 
We have also shown the persistence of spatiotemporal frag- 
mentation with respect to the zero order angular velocity 
(by considering both the MDI and the GONG data sets) , 
as well as under changes in the form of the factor f(9) 
that prescribes the latitudinal dependence of a (by con- 
sidering f(9) = sin 4 9 cos 9). We further found our model 
to be capable of producing butterfly diagrams which are in 
qualitative agreement with the observations. In this way 
we have found evidence that spatiotemporal fragmenta- 
tion is not confined to our original model only and that it 
can occur in more general dynamo models. 

Concerning our model, we should note that all our cal- 
culations were done with semi-open outer boundary con- 
ditions on the toroidal magnetic field (i.e. dB/dr = 0). 
In spite of relatively extensive searches, we have so far 
been unable to find spatiotemporal fragmentation with 
pure vacuum boundary conditions (B = 0). It is interest- 
ing to note in this connection that some examples of spa- 
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tiotemporal bifurcations in the literature (albeit in cou- 
pled maps) also have open boundary conditions (see e.g. 
Willeboordse & Kaneko 1994, Frankel et al. 1994). This 
may therefore be taken as some tentative support for the 
idea that spatiotemporal fragmentation may require some 
sort of open boundary conditions, or that at least it is 
easier to occur in such settings. 

We also made an extensive study of the spatial mag- 
netic field structure as well as the nature of dynamical 
variations in the differential rotation, including ampli- 
tudes and phases, as a function of depth and latitude. Our 
results demonstrate the presence of three main qualitative 
spatiotemporal regimes: (i) regimes where there is no de- 
formation of the torsional oscillations bands in the (r, t) 
plane, and hence no changes in phase or period through 
the convection zone; (ii) regimes where there is spatiotem- 
poral deformation in the oscillatory bands in the (r, t) 
plane, resulting in changes in the phase of the oscillations, 
but no changes in their period, and (hi) regimes with spa- 
tiotemporal fragmentation, resulting in changes both in 
phase as well as period/behaviour of oscillations, including 
regimes that are markedly different from those observed 
at the top, having either significantly reduced periods or 
non-periodic modes of behaviour. In all three cases, we 
found that torsional oscillations, resembling those found 
near marginal excitation, persisted above the fragmenta- 
tion level. 

These modes of behaviour can in principle explain 
a number of features that have been observed recently. 
These include: 

(a) Deformations of the type (ii) with or without frag- 
mentation of type (iii), can account for the reversal of 
the phase of the oscillatory behaviour above and below 
the tachocline, as reported by Howe et al. (2000a). 

(b) Spatiotemporal fragmentation can explain latitudi- 
nal dependence, as fragmentation occurs in both ra- 
dius/time as well as in latitude/time plots. In partic- 
ular, it can explain the possibility of oppositely signed 
tachocline shear at low and high latitudes, as has 
been found in observations by Howe et al. (2000a). 
To demonstrate this, we have plotted in Figs. ^9| and 

, the residuals of the differential rotation rate with 
respect to the radius and the latitude. 

(c) We observe the penetration of coherent torsional os- 
cillations into all the regions above the fragmentation 
level, with the penetration extending to slightly greater 
depths in the case of the GONG data. 

(d) Spatiotemporal fragmentation/bifurcation can in 
principle explain, purely dynamically, the relationship 
between the 11 year and possible 1.3 year oscillations 
near the top and bottom of the convection zone. In 
this connection, note that 11 = 2 3 x 1.3 years (com- 
patible with the result of Howe et al. 2000a being 
connected with three period halvings)! It can also ex- 




Radius 

Fig. 29. Variation of the perturbation to the zero order 
rotation rate in latitude and time, revealing the migrat- 
ing banded zonal flows, after the transients have died out, 
using the MDI data. Observe that residuals with different 
signs can occur over small latitude or radii bands. Param- 
eter values are as in Fig. 
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Fig. 30. Variation of the perturbation to the zero order 
rotation rate in latitude and time, revealing the migrating 
banded zonal flows, after the transients have died out. 
Parameter values are as in Fig. ||. 



plain non-periodic dynamical behaviour (compatible 
with the findings of Antia & Basu 2000). 

Finally, apart from providing a possible theoretical 
framework for understanding such phenomena, this sce- 
nario could, by demonstrating the different qualitative dy- 
namical regimes that can occur in the dynamo models, 
also be of help in devising strategies for future observa- 
tions. 
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